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Abstract 

Growth (and resorption) of biological tissue is formulated in the continuum setting. The treatment 
is macroscopic, rather than cellular or sub-cellular. Certain assumptions that are central to classical 
continuum mechanics are revisited, the theory is reformulated, and consequences for balance laws and 
constitutive relations are deduced. The treatment incorporates multiple species. Sources and fluxes 
of mass, and terms for momentum and energy transfer between species are introduced to enhance the 
classical balance laws. The transported species include: (i) a fluid phase, and (ii) the precursors and 
byproducts of the reactions that create and break down tissue. A notable feature is that the full extent 
of coupling between mass transport and mechanics emerges from the thermodynamics. Contributions 
to fluxes from the concentration gradient, chemical potential gradient, stress gradient, body force and 
inertia have not emerged in a unified fashion from previous formulations of the problem. The present work 
demonstrates these effects via a physically-consistent treatment. The presence of multiple, interacting 
species requires that the formulation be consistent with mixture theory. This requirement has far-reaching 
implications. A preliminary numerical example is included to demonstrate some aspects of the coupled 
formulation. 



1 Background 

Development of biological tissue, when described in the biomechanics literature, is generally broken down into 
the distinct processes of growth, remodelling and morphogenesis. Growth, or conversely, resorption, involves 
the addition or loss of mass. Remodelling results from a change in microstructure, which could manifest 
itself as an evolution of macroscopic quantities such as state of internal stress, stiffness or material symmetry. 
It also appears sometimes as fibrosis or hypertrophy. Morphogenesis involves both growth and remodelling, 
as well as more complex changes in form. A classical example of morphogenesis is the developm ent of an 



embryo from a fertilized egg. These terms are based on the definitions developed by Taber (1995), and will 
be followed in this work. In the present communication we will focus exclusively on growth, its continuum 
formulation, and the implications that the process h olds for the standard ma chinery of continuum mechanics. 



Remodelling is treated in an accompanying paper I Garikipati et al. . 2003|) . The larger and more complex 



problem of morphogenesis will not be treated in this body of work. 
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The ideas here are applicable to soft (e.g., tendon, muscle) and hard (e.g., bone) tissue. In this paper, 
growth of biological tissue will be treated at a macroscopic scale. The continuum formulation (e.g., consti- 
tutive laws) at this scale may be motivated by cellular, sub-cellular or molecular processes. However, we 
will not explicitly model processes at this fine a scale. The formulation can be applied with a specific tissue 
e.g., muscle, as the body of interest. Our experim ents, described se parately, are on self-organizing tendon 



i Calve et al 



and cardiac muscle constructs IjBaar et al 



engineered in vitro. 



The principal notion to be borne in mind while developing a continuum formulation for growth is that one 
is presented with a system that is open with respect to mass. Scalar mass sources and sinks, and vectorial 
mass fluxes must be consid ered. A mass source was first introduced in t he context of biological gro wth by 
Cowin and Hegediisl 1 197fi[) . The mass flux is a more recent addition of Epstein a nd Maugftl iboOf j ) who 
however, did not elaborate on the specific nature of the transported species. iKuhl and Steinma nn (2002) 
also incorporated the mass flux and specified a Fickean diffusive constitutive law for it. In their paper the 
diffusing species is the material of the tissue itself. The approach to mass transport that is followed in our 
paper is outlined in the next two paragraphs. 

In order to be precise about the physiological relevance of our formulation, we have found it appropriate 
to adopt a different approach from the papers in the preceding paragraph in regard to mass transport. We 
do not consider mass transport of the material making up any tissue during growth. Instead, it is the 
nutrients, enzymes and amino acids necess ary for growth of tissue, byproducts of metabolic reactions, and 
the tissue's fluid phase i Swartz et all 1999h that undergo diffusion 1 in our treatment. There do exist certain 
physiological processes in which cells or the surrounding matrix migrate within a tissue. One such process 
is observed when leukocytes (white blood cells) such as neutrophils and monocytes are signalled to pass 
through a capillary wall and are in duced, by specific chemical attractors, to migra te to a site of infection. 
This is the process of chemotaxis (jGuvton and Hail ll 99(1 Iwidmaier et all hoO.*^ . The migrant cells or 
matrix then participate in some form of cell proliferation or death. Fibroblasts also migrate within the extra 
cellular matrix during wound healing. A third example is the migration of stem cells to different locations 
during the embryonic development of an animal. These processes involve very short range diffusion, and can 
be treated by the appr oach described in this paper. We have chosen to focus upon homeostasis, defined by 
Widmaier et al. 1 200.4 ) as "... a state of reasonably stable balance between the physiological variables. . . " 2 . 
Since, to the best of our knowledge, processes of the type just described are not observed during homeostatic 
tissue growth, we will ignore transport of the solid phase of the tissue. 

The processes of cell proliferation and death, hypertrophy and atrophy, are complex and involve several 
cascades of biochemical reactions. We will treat them in an elementary fashion, using source/sink terms that 
govern inter-conversion of species, and the mass fluxes that supply reactants and remove byproducts. The 
treatment will be mathematical; specific constituents will not be identified with any greater detail than to say 
that they are either the tissue's solid phase, the interstitial fluid phase, precursors of the solid phase (these 
would include amino acids, nutrients and enzymes), or byproducts of reactions. We will return to explicitly 
incorporate biochemical and cellular processes within our description of mass transport in a subsequent 
paper. 

Vir tually all biological tissue consists of a solid and fluid phase and can be treated in the con text of mixture 
theory ( Truesdell and TouoirJ . 196ot Truesdell and Noll 196.4 Bedford and Drumheller , 198.4) . When growth 



iWe use the terms "ma ss transport" and "diffusion" interchangeably. 

i Widmaier et alJ 1 2003) go on to say: "This simple definition cannot give a complete appreciation of what homeostasis truly 
entails, however. There probably is no such thing as a physiological variable that is constant over long periods of time. In fact, 
some variables undergo fairly dramatic swings about an average value during the course of a day, yet may still be considered 
'in balance'. That is because homeostasis is a dynamic process, not a static one." 
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is of interest, additional species (reactants and byproducts) must be considered as outlined above. The solid 
phase is an anisotropic composite that is inhomogeneous at microscopic and macroscopic scales. The fluid, 
being mostly water, may be modelled as incompressible, or compressible with a very large bulk modulus. 
This level of complexity will be maintained throughout our treatment. The use of mixture theory leads to 
difficultie s associated with part itioning the boundary traction into portions corresponding to each species. 



Raiagooal and Winemar] l)199ff ). suggested a resolution to this problem that holds in the case of saturated 
media — a condition that is applicable to soft biological tissue. An alternativ e is to app l y the theory of porou s 



media that grew out of the classical work of Fick and Darcy in the 1800s jTerzaehil Il943t 



de Boei 



ijorous 
2OO0l) . 



In this approach fluxes are introduced for each species. Since a species that diffuses must do so within some 
medium, one may think of the various constituents diffusing through the solid phase 3 . This strategy has 
been adopted in the present work. 



1.1 Recent work 



In a s implification that avoided the complexity of mixture theory or porous media, Cowin and Heeedusl 
(|l9T6^) accounted for the fluid pha se via irreversible sou r ces an d fl uxes of momentum, energy a nd entropy. 
This approach was also followed by Epstein and Maueinl ( 200o|) and 



Kuhl and Stcinman: 



= 



2002(1 . While the 

approach followed in the present paper, i.e. derivation of a mass balance law with mass source and flux, 
and postulating sources an d flu xes for momentum, energy and entropy, has been attempted recently by 



Epst ein and M augii 



ul ll200(T|) . 



and 



Kuhl and Steinmanr 



there are important differences between those 



works and our paper. Epstein and Maugin conclude that the mass flux vanishes unless the internal energy 
depends upon strain gradient terms (a second-order theory). This view ignores Fickean diffusion (where 
the flux is linearly dependent upon concentration of the relevant species). Our treatment also results in the 
dependence of mass flux upon strain gradient, but without the requirement of a strain gradient dependence 
of the internal energy. Instead, the dissipation inequality motivates a constitutive relation for the mass flux 
of each transported species. When properly formulated in a thermodynamic setting, the mass flux can be 
constrained to depend upon the strain/stress gradient and the gradient in concentration (mass per unit 
system volume) of the corresponding specie s. The latter term is the Fic kean contribution to the mass flux. 
The form obtained is essentially identical to lde Groot and Mazurl 1,1984). 



While modelling hard biological tiss ue, it is common to assume a Fickean flux, and a mass source that 



depends upon the strain energy density (Harrigan and Hamilton! . 119 93). and therefore upon the strain. This 
introduces coupling between mass transport and mechanics. One possible difficulty with this approach is 
that one could conceive of a mass source that satisfies other requirements, such as the dissipation inequality, 
but does not depend upon any mechanical quantities. Strain- or stress-mediated mass transport would then 
be absent in the boundary value problems solved with such a formulation. 

The present paper is aimed at a complete treatment of mass transport, coupled with mechanics, for the 
growth problem. Initial sections (Sections EHEJ) treat the balance of mass, balance of linear and angular 
momenta, the forms of the First and Second Laws for this problem, and kinematics of growth, respectively. 
The Clausius-Duhem inequality and its implications for constitutive relations are the subject of Section H3 
Examples are provided as appropriate to illustrate the important results. A preliminary numerical example 
appears in Section A discussion and conclusion are provided in Section |H1 

3 A more sophisticated, and physiologically-correct, description is that the interstitial fluid diffuses relative to the solid phase, 
while precursors and byproducts of reactants diffuse with respect to the fluid. 
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2 Balance of mass for an open system 



The body of interest, 23, occupies the open region fio C R 3 in the reference configuration. Points in 23 are 
parameterized by their reference positions, X. The deformation of 23 is a point-to-point map, tp(X,t) £ R 3 , 
of £!o, carrying the point at X to its current position x = ip(X,t), at time t £ [0, T]. In its current 
configuration, 23 occupies the open region f2 t = cp t (fl(y), fl t C K 3 (see Figure^. The tangent map of tp is 
the deformation gradient F := d(p/dX. 




Figure 1: The continuum body with diffusing species and under stress. 



The body consists of several species, of which the solid phase of the tissue is denoted s, and the fluid 
phase is f. The remaining species, a, . .. ,u>, are precursors of the tissue or byproducts of its breakdown 
in chemical reactions. The index i will be used to indicate an arbitrary species. (Where appropriate, we 
will use the term "system" to refer to the species collectively. Where the presence of several species is 
an unimportant fact, we will use the term "body". ) The system is open with respect to mass. Species 
a, . . . , u> have sources/sinks, IT", . . . , IF', and mass fluxes, M a , . . . , , respectively. The sources specify 
mass production rates per unit volume of the body in its reference configuration, Qq. The fluxes specify 
mass flow rates per unit cross-sectional area in fV Importantly, these flux vectors are defined relative to the 
solid phase. As discussed in Section^ the solid tissue phase has only a mass source/sink associated with it, 
and no flux. The fluid tissue phase has only a mass flux, and no source/sink 4 . 

Since the solid phase, s, does not undergo transport, its motion is specified entirely by ip(X,t). We 
describe the remaining species f, a, . . . , u> as convecting with the solid phase and diffusing with respect to it. 
They therefore have a velocity relative to s. Since the remaining species convect with s, it implies a local 
homogcnization of deformation. The modelling assumption is made that at each point, X, the individual 
phases undergo the same deformation. 

We define concentrations of the species p L = p L f L as masses per unit volume in f^o- The intrinsic species 

density is p L , and f L is the volume fraction of l, for t = s, f, a, . . . , u>. In an experiment it is far easier to 

measure the concentration, p L , rather than the intrinsic species density, Pq 5 . The concentrations also have the 

property ^2 Po — Pcb the total material density of the tissue, with the sum being over all species s, f, a, . . . , u>. 
l 

4 Considering the case of the lymphatic fluid, this implies that lymph glands are assumed not to be present. 
5 In our experiments we have m easured the mass c oncentration p L Q of collagen in engineered tendons grown in vitro. These 
results will be presented elsewhere <Calve et allboO-A . 
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The concentrations, p l , change as a result of mass transport and inter-conversion of species, implying that 
the total density in the reference configuration, po, changes with time. They are parameterized as Pq(X , t). 



2.1 Balance of mass in the reference configuration 

Recall that flo is a fixed volume. The statement of balance of mass of the solid phase of the tissue, written 
in integral form over Oq, is 



^ J P l(X,t)dV = J IL s (X,t)dV. 



(1) 



Since the solid phase of the tissue does not undergo mass transport, there is no associated flux. Localizing 
the result gives 



8t^ V 
where the explicit dependence upon position and time has been suppressed. 

The fluid phase of the tissue, f , may be thought of as the interstitial or lymphatic fluid that perfuses the 
tissue. As explained above, we do not consider sources of fluid in the region of interest. The fluid therefore 
enters and leaves fio as a flux, M f . The balance of mass in integral form is 

£ / pl(X, t)dV = - f M f (X, t) ■ NdA, (3) 

n 

where N is the unit outward normal to the boundary, dflo- Applying the Divergence Theorem to the surface 
integral and localizing the result gives 

d -§ = -V • M f , (4) 

where V(») is the gradient operator defined on Sloj and V • (•) denotes the divergence of a vector or tensor 
argument on Qq. 

For the precursor and byproduct species, l = a, . . . ,u>, the balance of mass in integral form is 

^ J P L (X, t)dV = J U L (X, t)dV - J M\X, t) ■ NdA. (5) 

In local form it is 

M =n i - V-M l , Vi = a,...,u. (6) 

Of course, this last equation is the general form of mass balance for any species t, recalling that in particular, 
M s = and II = 0. This form will be used in the development that follows. 

The fluxes, M L , Vt = f, a, . . . ,ui, represent mass transport of the fluid, of precursors to the reaction 
site, and of byproducts from sites of tissue breakdown. The sources, IT , Vt = s, a, . . . , arise from inter- 
conversion of species. The sources/sinks in © and © are therefore related, as tissue and byproducts are 
formed by consuming precursors (amino acids and nutrients, for instance). To maintain a degree of simplicity 
in this initial exposition, we will restrict our description of tissue breakdown to the reverse of this reaction. 
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The sources, IP for various species, satisfy a relation that is arrived as follows: Summing 10 over all 
species leads to the law of mass balance for the system, 



^ dt 



-J^Y^IF-V-M 1 ), (7) 



where the sum runs over all species, with n f = and M b = 0. Alternatively, in writing the mass balance 
equation for the system, the interconversion terms (sources/sinks) play no role, and only the fluxes at the 
boundaries need be accounted for. In integral form we have 



Applying the Divergence Theorem and localizing leads to 

£^ = £(-™')- (8) 

Comparing the equivalent forms Q and JBJ it emerges that the sources and sinks satisfy 

vjir = o, (9) 



a conclusion that is consistent with classical mixture theory ( Truesdell and Noll 1965^) . The section that 



follows contains an example in which the Law of Mass Action is invoked to describe a set of inter-related 
sources, 11'. 



2.1.1 Sources, sinks and stoichiometry: An example based upon the Law of Mass Action 

The conversion of precursors to tissue and the reverse process of its breakdown are governed by a series of 
chemical reactions. The stoichiometry of these reactions varies in a limited range. Continuing in the simple 
vein adopted above, it is assumed that the formation of tissue and byproducts from precursors, and the 
breakdown of tissue, are governed by the forward and reverse directions of a single reaction: 

E^W-N- (io) 

i— a 



Here, n L is the (possibly fractional) number of moles of species t in the reaction. For a tissue precursor, 
n L > 0, and for a byproduct, n L < 0. By the Law of Mass Action for this reaction, the rate of the forward 
reaction (number of moles of s produced per unit time, per unit volume in Qq) is kf Y[ [Po]™'; where J~] 

i— a 

on the right hand-side denotes a product, not to be confused with the source, II. The rate of the reverse 
reaction (number of moles of s consumed per unit time, per unit volume in Q ) is & r [PoL where kf and k r 
are the corresponding reaction rates. Assuming, for the purpose of this example, that the solid phase is a 
single compound, let the molecular weight of s be M s . From the above arguments the source term for s is 

n s - Un^r-fcr^])M s , (ii) 



G 



Since the formation of one mole of s requires consumption of n L moles of t, we have 

ir = -(k f[[p^ -fc r [s]| n t M t , (12) 



0=a 



where M t is the molecular weight of species i. Since, due to conservation of mass, M s = ^'^VT^ the sources 

l— a 

satisfy £ W = 0. 



l— s,a 



2.2 Balance of mass in the current configuration 

In the current configuration, fi t , the concentration, source and mass flux of species l are p L {x,t) 1 7r l (a;,i), 
and m L (x,t) respectively. The boundary is <9f2 t , and has outward normal n (Figure Since the de- 
formation, ip(X,t), is applied to the body (the system), standard arguments yield p l = ^(detF)" 1 , and 
tt l = ir^det-F 1 ) -1 . By Nanson's formula, the normals satisfy nda — F ~ T NdA, where dA and da are the area 
elements on the boundaries <9f2o and dttt respectively. The Piola transform then gives m L = (detF)~ 1 FM l . 

The balance of mass in integral form for the solid tissue phase is 

d 
dt 



P S {x,t)dv = j TT S {x,t)dv. (13) 

u t n t 



Applying Reynolds' Transport Theorem to the left hand-side, localizing the result and employing the product 
rule gives, 

V xP s +p s V x -v = ir s , (14) 

at 

or, 

^ - 7T S - p s V x ■ v, (15) 

where v(x,t) = V(X, t) oif^ 1 {x, t) is the spatial velocity, V x («) is the spatial gradient operator, and V x • (•) 
is the spatial divergence. The time derivative on the left hand-side in (|15|l is the material time derivative, 
that may be written explicitly as (dp s (x,t)/dt)x, implying that the reference position is held fixed. 



For the fluid tissue phase, f , the integral form is, 

■ J p f (x,t)dv = — J m i (x,t)-nda. (16) 



dl 



n t do. t 

Invoking the Divergence Theorem in addition to the arguments used above for the solid phase now gives 

^ = -V x ■ m f - p l V x ■ v. (17) 

Finally, for the precursor and byproduct species, t = a, ...,u>, proceeding as above gives the local form 
in Q t : 

% = 7T l - V x ■ m l - p'V x • v, Vi = a,...,w. (18) 
dt 
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3 Balance of linear and angular momenta 



3.1 Balance of linear momentum 

We first write the balance of linear momentum in the reference configuration, Q,q. The body, S, is subject 
to surface traction, T, and body force per unit mass, g. The natural boundary condition then implies that 
T = J2 L P^N on dfto, where P L is the partial first Piola-Kirchhoff stress tensor corresponding to species 
i, and the index runs over all species. Thus, P l N is the corresponding partial traction. The mass fluxes, 
M L , (t = f, a, . . . , cj), and mass sources, IT , (t = s, a, . . . , uS) make important contributions to the balance 
of linear momentum, as shown below. 

The body undergoes deformation, (p(X,t), and has a material velocity field V(X,t) = dip(X,t)/dt. In 
discussing momentum and energy, it proves convenient to define a material velocity of species t relative to 
the solid phase as V L = (1/ p L Q )FM L . Recall (from Section |2J) that the remaining species are described as 
deforming with the solid phase and diffusing relative to it. Therefore F is common to all species. The spatial 
velocity corresponding to V 1 is v l — (1/ p L )m L = V L , by the Piola transform. Since fluxes are defined relative 
to the solid tissue phase, which does not diffuse, the total material velocity of the solid phase is V, and for 
each of the remaining species it is V + V L , t, = f, a, . . . , u). Formally, we can write the material velocity as 
V + V\ L = s,f, a, . . . , u> with the understanding that V s = 0. Likewise, II f = 0. This convention has been 
adopted in the remainder of the paper. 

The balance of linear momentum of species i written in integral form over Qq is, 

^ J ftty + v^dv = J P L gdV + J P L Q q l dv + 1 n'(v + v')dv 



J P L NdA - J (V + V L )M L ■ NdA, (19) 



8Qq 8Qq 



where g is the body force per unit mass, and q L is the force per unit mass exerted upon t by the other species 
present. Attention is drawn to the fact that the mass source distributed through the volume, and the influx 
over the boundary affect the rate of change of momentum in (|19|) . Summing over all species, the balance of 
linear momentum for the system is obtained: 



E 



^ / p L (V + V^)dV = E/ Po9W + Y^ J Pol^V 

u ^ v + vL ) dv +E / pLNdA 



1 dCl 



J (V + V L )M L ■ NdA, (20) 



1 an 

The interaction forces, PoQ 1 , satisfy a relation with the mass sources, IT', that is elucidated by the 
following argument: The rate of change of momentum of the entire system is affected by external agents 
only, and is independent of internal interactions of any nature (q L and IT). This observation leads to the 
following equivalent expression for the rate of change of linear momentum of the system: 
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Y,^ t fpo( V + V ') dV = Jp°9dV+ J PNdA 

-^2 J (V + V L )M' ■ NdA. (21) 

Here, P = ^ P L and po ~z~2Po- Since both and (|21[) represent the balance of linear momentum of the 
system, it follows that, 

£ J pW&v + ]T | n 4 (v + v 1 )^ - o (22) 

Recalling the relation between the sources and localizing leads to 



£( P ^ + irn = o, 



(23) 



a result that is also consistent with classical mixture theory I Truesdell and Noll 1965^ . 



Having established we return to the balance of linear momentum for a single species (|19f) in order 
to simplify it. Writing (V + V')[M L ■ TV] as ((V + V) <g> M L )N , and using the Divergence Theorem, 



+ J (n l (v + v L ) + v • p l ) dv 

-J V -((V + V L )®M L )dV 



(24) 



Using the mass balance equation Jfljl, and applying the product rule to the last term gives 
J p' ^ (V + V) dV = Jti(g + q L ) dV 

+ J (V • P L - (V (V + V")) AT) dV 
n 



(25) 



Localizing this result gives the balance of linear momentum for a single species in the reference configuration: 



P '°di {V + =p ^ 9 + ql) + V • P' - (V (F + V')) M< 



(26) 
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The balance of linear momentum for a single species in the current configuration, O t , is obtained via 
similar arguments and the Reynolds Transport Theorem: 

d 

p b —{v + v L ) = p L {g + q L ) + V x -(T L 

- (V x (v + O) m< - / (V a (v + v 1 )) v, (27) 
where er L — (detF)~ 1 P l F T is the partial Cauchy stress of species i. 

3.2 Angular Momentum 

For the purely mechanical theory, balance of angular momentum implies that the Cauchy stress is symmetric: 

T 

<j L = <j l . We now re-examine this result in the presence of mass transport. At the outset one might expect 
a statement on the symmetry of some stress or stress-like quantity. We derive this result for any species, t, 
beginning with the integral form of balance of angular momentum written over Oo. 

^J<pxrf>(V + V l )dV = Jtpx [p&( ff + g 4 )+n 4 (V + V l )]dV r 

+ J <p x {P L - (V + V L ) ®M L )NdA (28) 
Applying properties of the cross product, the Divergence Theorem and product rule gives 

J vxp^ + cpx f^(y + v l ) + ^{y + v l )jdv = 

J<pxp L Q (g + q L + U L (V + V"))dV 

o„ 

+ J {<p x V • P L - ip x (V (V + V L ) M L )) dV 

n 

J (-if x (V + V L ) V ■ (M 1 )) dV 
o„ 

- J e: ((P l - (V + V L ) ® M L ) F T ) dV, (29) 
o 

where e is the permutation symbol, and e: A is written as Sijk-Ajk in indicial form, for any second-order 
tensor A. Using the mass balance equation |(JjJ|, and balance of linear momentum (I2fi|) . we have 

J V x P ' V L dV = - J e: \p l - ( V + V ) <8> F T dV. 

o oo \\ p^F-'Vy J 



Oo 
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Recalling the relation of the permutation symbol to the cross product, and the indicated relation between 
M L and V L leads to 

= - J e: ((P l - V 1 <& p^F^V") F T ) dV. (30) 

n 

Localizing this result and again applying the properties of the permutation symbol we are led to the symmetry 
condition, 

{P b -V L ® p'qF^V) F t = F (P l -V b ® ^F _1 y t ) T . (31) 

But, (V <g> F~ 1 V-)F T = V L (g> V\ Thus, the symmetry P l F T = F(P l ) T that results from conservation of 
angular momentum for the purely mechanical theory, is retained in this case. The partial Cauchy stresses 

T 

are therefore symmetric: <t' = cr l . This is in contrast with the non-symmetric Cauchy stress arrived at by 
Epstein and Mauginl ( 200o|) . The origin of this difference lies in the fact that these authors use a single species 



with V — d<p/dt as the material velocity, rather than multiple species with material velocities V + V L . 

4 Balance of energy and the entropy inequality 
4.1 Balance of energy 

Since mass is undergoing transport with respect to 23, and inter-conversion between species i = s, a, . . . , w, 
it is appropriate to work with energy and energy-like quantities per unit mass. In addition to the terms 
introduced in previous sections, the internal energy per unit mass of species i is denoted e' ; the heat supply 
to species i per unit mass of that species is r l ; and the partial heat flux vector of i is Q L , defined on ^o- An 
interaction energy appears between species: The energy transferred to i by all other species is e L , per unit 
mass of l. In the arguments to follow in this section we will use the fluxes M L and the associated velocities, 
V L . Working in Oo, we relate the rate of change of internal and kinetic energies of species I to the work done 
on l by mechanical loads, processes of mass production and transport, heating and energy transfer: 

^/po U + \\\V + VT) dV = J (ffa-(V + V) + tir>)dV 

+ J pW ■ (v + v L )dv 
+ J fir ^ + I||v + vi a )+^ dv 



J ({V + V u ) P L M L (e l + i||V + V L f\j - Q h \ ■ NdA. 



(32) 



ao 

Summing over all species, the rate of change of energy of the system is, 



11 



(e L + \\\V + V<-f) dV = Y,J Ote-(V + V) + tir>)dV 

+ E/ • (V + V')dV 

+ E / (n l (V + i||v + vi 2 )+^ dF 

2 / f ( V + V 1 ) • P l - M' Ye* + i||V + V 4 || 3 ^ -Q^j -NdA. 



(33) 



The inter-species energy transfers are related to interaction forces and mass sources. To demonstrate this, 
we proceed as follows: The rate of change of energy of the system can also be expressed by considering the 
system interacting with its environment, in which case the internal interactions between species (interaction 
forces, mass interconversion and inter-species energy transfers) play no role. This viewpoint gives, 

E^fo (e l + l\\V + VT) dy = E/ (Po9-(V + V l )+p^)dV 
+ E / {{V + V L )-P L -M L [e' + \\\y + V 4 || 3 ^ -Q^ -NdA. 



(34) 



ao 

Since l|3"3j) and (|3"4")l are equivalent, it follows that, 

E \J^PW-(V + V')+W [e" + 1 -\\V + V' 
and on localizing this result, 



L II 2 \ I _£ ~i 



p^ )dV = 0, (35) 



^^.(v + v^ + w^ + ^wv + vt) ) =')• ^>> 

This result relating the interaction energies to interaction forces betwe en species, their sources a nd relative 
velocities, is identical to that obtained from classical mixture theory i Truesdell and Noll 1965b . Together 



with 10 and (|23|) it demonstrates that the present formulation is consistent with mixture theory. 

Equation (|32[) for the rate of change of energy of a single species can be further simplified by applying 
the Divergence Theorem and product rule, giving first, 
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/ (M (e- + \\\V + V|") + 4 (e' + + V|")) & = 

/ (po9 • (V + V 1 ) + p&r* + IT (V + i|| V + + pfce 4 ) dU 

+ J p' q L -(V + V L )dV 
+ J ((V + V L ) ■ V • P L + P' : V (V + V 1 )) dU 

-| (v-(M') ^ + i||y + v*f^dV 

n 

- y ((Ve l + (V + V ) • V {V + V L )) ■ (M l ) - V ■ Q L ) dV. (37) 

Using the balance of mass ©, balance of linear momentum (|26f) . and localizing the result, we have, 
de L 

p& — = P L :V(V + V L )-V-Q'+ P y+p' e'-Ve'-M" (38) 
Summing over i gives, 

t 

(P K - F + P U VV - V • Q l + p L a r L + p L l L - Ve' ■ M') (39) 

L 

Substituting for ^ p L Q e L from <j5B)t . 

E' 9 ^ = ^ (P'iF + P 1 : W l - V-Q'+p^r 1 - Ve l • M l ) 

- £ • (v + v L ) n- (V + \\\v + vi 2 ) ) . (40) 

This form of the balance of energy is most convenient for combining with the entropy inequality leading 
to the Clausius-Duhcm form of the dissipation inequality. 
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4.2 The entropy inequality: Clausius-Duhem form 



Let rf be the entropy per unit mass of species i, and 9 the absolute temperature. The entropy production 
inequality holds for the system as a whole. Accordingly, we write 



Po 7 



dV 



J yM L ■ Nrf + Q--N^ dA. 



(41) 



1 aa 



Applying the Divergence Theorem, using the mass balance equation ©, and localizing the result, we have 
the entropy inequality, 



NT ' ®^ ~> ST (I* 1 * T7 » ^ v 91 j. ve ■ 91 



(42) 



Now, multiplying Equation (|42|l by 9, subtracting it from Equation l|4(J|l and using l|2tj[l for p^q 1 gives, 



+ Y,(p l o§- t (V + V')-pl l g-V-P L + V(V + V')M^-(V + V L ) 
- Yl (P u - F - P L VV + (Ve l - 9Vri l ) ■ M'j < 



(43) 



Equation l|43[) is the reduced entropy inequality — also referred to as the Clausius-Duhem inequality — for 
growth processes. 



5 The kinematics of growth 

The formulation up to this point has introduced some elements of coupling between mass transport, me- 
chanics and thermodynamics. Mass transport and mechanics are further coupled due to the kinematics of 
growth. Local volumetric changes take place as species concentrations evolve. As concentration increases, 
the material of a species swells, and conversely, shrinks as concentratio n decreases. This observation has 
led to an active field of study within the literature on biological growth Jskalak L 



Taber and Humphrey . 2001 ; Lubarda and Hoger , 20021: 



Ambrosi and Mollic 



3EI 



1981 



Skalaketal 



Our treatment follows 



in the same vein. 
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5.1 The elasto-growth decomposition 



Finite strain kinematics treats the total deformation gradient as arising from a geometrically-necessary 
elastic deformation accompanying growth, as well as a separate elastic deformation due to an external stress. 
The defor mation grad i ent is subject to a split reminiscent of the classical decomposition of multiplicative 



plasticity l|Bilbv et al 



1957 



subject to 
Leelll969h 



At a continuum point the reference concentration of each species admits the notion of an "original" state 
in which the concentration of a species is p L mg (X). This is a state that may never be attained in a physical 
system. However, if attained, the corresponding species would be stress-free in the absence of deformation. 
Neglecting other possible kinematics (such as plasticity) and microstructural details, the set of quantities 
{Pq, . . . ,Pq}, and the temperature, 9, fully specify the reference state of the material at a point. As mass 
transport alters the reference density to its value p L {X,t), the species swells if p L Q > p L mg , and shrinks if 
Pa < Porg- Assuming that these volume changes are isotropic leads to the following growth kinematics: For 
each species, one can define a "growth deformation gradient tensor", F 6 := -^-1, where 1 is the second- 

Porg 

order isotropic tensor. The tensor F g is analogous to the plastic deformation gradient of multiplicative 
plasticity. As the ratio Po//°org i s a local quantity, F 6 varies pointwise and adjacent neighborhoods will, in 
general, be incompatible due to the action of F g alone. However, further elastic deformation, F occurs 
to ensure compatibility, leading to an internal stress, that can in general be different for each species. The 
action of these kinematic tangent maps can be conceived of in the absence of external stress. With an 
external stress, there is further elastic deformation, F° , common to all species. This sequence of maps is 
pictured in Figure [3 The kinematic relations are: 

F = F c F°'f s \ F s ' = Jp-1. (44) 

Porg 

Clearly, the elastic deformation gradients can be combined to write F e = F°F , the "total" elastic defor- 
mation gradient of species i. 



6 Restrictions on constitutive relations from the Clausius-Duhem 
inequality 

As is the practice in field theories of continuum physics, we use the Clausius-Duhem inequality l|43|) to obtain 
restrictions on the constitutive relations. We begin with very general assumptions on the specific internal 
energy of each species: e L — e'(_F c ,i] L ,p L ). Then, applying the chain rule and regrouping some terms, (|43() 
becomes 

*f \dF c ) £f ) at 



J2(p' F- T (Ve' -6V,f))-V 
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Figure 2: The kinematics of growth. 



£ir ( e i + i||v + v*| 



ve-Q L 



(45) 



Inequality i|43|) represents a fundamental restriction upon the physical processes during biological growth. 
Any c onstitutive relations that are prescribed must satisfy this restriction, as is well-known ( Truesdell and Noll 
196.51) . Guided by l|43|) . we prescribe the following constitutive relations in classical form: 



(46) 



de' 
dr\ L 



(47) 
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~ L 



{ pl °^r ~ pl ° 9 ~ v ' pL + (w) ML + P * F ~ T (Vet ~ 9Vlf) ' f - 1 4Si 

where w ■ Dw > 0, el 3 



Q l = -K L V6, w ■ K u w > 0, V«)€l 3 (49) 



With the constitutive relations H46H49fl ensuring the non-positiveness of certain terms the entropy inequality 
is further reduced to 

+E (po v ' ■ {^f + F ~ lvi ) + nt ( e ' + \\\y + y L \\ 2 )) 

+ E (pa^iV + V^-Pod- V-P" + V(V + V L )M^ ■V<0. (50) 
The left hand-side of H5()(l is the dissipation, D, a quantity that we will return to below. 

lT 

Equation (|46(1 specifies a constitutive relation for P L F & , which is a truly elastic stress. Equation (|47(1 
implies a uniform temperature in each species, and corresponds with the definition usually employed for 
temperature in thermal physics. The heat flux in species i is given by the product of a positive semi- 
definite conductivity tensor, K L and the temperature gradient, which is the Fourier Law of heat conduction. 
Equation Ij48(l requires more detailed discussion which appears below. 

6.1 Constitutive relations for fluxes 

Turning to l|48|) . we first point out that since AT = p'qF-W-, this is an implicit relation for V L . Rewriting 
it as an explicit one for PqV l we have, 

{, frVVF^Y 1 D L 

p ° v = + 7* 

fii^- - A/ - V • P + p'„F- T (Ve' - SV,f )) | . (51) 
The constitutive relation for flux, M l = p L F~ 1 V is then obtained: 

i / b L WF~ 1 \ 1 b L T 

M =F- 1 \1 + ■ —F~ T 



Po J Po 
D 
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PoF T — - p' F T g - F T V ■ P l + & (Ve 4 - 0Vt, 1 ) 



> . 



(52) 



As is common in descriptions of mass transport, the tensor delineated as D L will be referred to as the 
mobility tensor of species i. Recall that in (|48|l we have taken D to be positive semi-definite 6 . The flux 
is thus written as the product of a mobility tensor and a thermodynamic driving force, CP. This is the 
Nernst-Einstein relation. We proceed now to examine the four separate terms in the thermodynamic driving 
force: 

T = -P U F T — + P L Q F T g + F T V ■ P L - p L Q (Ve l - OV^) (53) 

The first two terms respectively represent the influences of inertia and body force. Thus, the inertial effect is 
to drive species i in the opposite direction to the body's acceleration. The body force's influence is directed 
along itself. The third term represents the stress divergence effect. In the case of a non- uniform partial 
stress, P L , there exists a thermodynamic driving force for transport along P\ We demonstrate this effect 
for the case of the fluid species in Section RT2l below, for which it translates to the more intuitive notion of 
transport along a fluid pressure gradient. 

The fourth term in CP admits the following interpretation: The Legendre transformation ip l = e L — 6rf 
allows one to rewrite Ve' — QWr) L as 'Vip'\s (at uniform temperature), where ip L is the mass-specific Hclmholtz 
free energy. An assumption inherent in the development that began in Section [5] is that any mass entering 
or leaving f2o at a point X on the boundary, dilo, has the field values p L ,e L ,r]'', 8, and ip L corresponding to 
X . Likewise, the incremental mass of species t created or absorbed via the source/sink IP at X has the 
field values of that point. Consider a sufficiently small neighborhood of a point, say N(X) C flo- Changing 
the mass of species i in N(X) by i5m' units causes a change in the Helmholtz free energy of i in N(X) by 
5^! L = ip L 5m 1 '. By definition therefore, if) 1 — d^ L /dm 1 . This derivative gives the chemical potential, p l , of the 
transported species, l. Thus, we have p L = e' — drf , and Ve l — 9"Vf]' = V/i'|e. This last term in CP thus 
represents the thermodynamic driving force due to a chemical potential gradient. 

It has rec ently come to our attention that the constitutive relation for flux (|52|l is precisely the result 



arrived at bv Ide Croot and Mazurl (|1984ri . including the identification of the chemical potential gradient 
term. However, their approach involves a slightly different application of the Second Law, and a less detailed 
treatment of the mechanics. 

The gradient of internal energy in (|53(l leads to a strain gradient-dependent term. A concentration 
gradient-driven term arises from the gradient of mixing entropy. Together with the other terms that were 
remarked upon above, they represent a complete thermodynamic formulation of coupled mass transport and 
mechanics. This is the central result of our paper. 

6.2 Transport of the fluid species: The example of an ideal fluid 

Consider the stress divergence term i r ' T V • P\ An elementary calculation gives 

F T V P' =V (^F T P^ - VF T : P\ (54) 
6 If IIZJVVF-VpoII << 1, we have D L £>'. 
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In indicial form, where lower/upper case indices are for components of quantities in the current/reference 
configuration respectively, this relation is 

F iKPiJ,J = ( F iKPij)j ~ FiKjPiJ- 

For an ideal fluid, supporting only an isotropic Cauchy stress, pl, we have P l = dct(F)pF~ T , where p is 
positive in tension. The arguments that follow assume this case. (The more general case of a non-ideal, 
viscous fluid will merely have additional terms from the viscous Cauchy stress.) The stress divergence term 
is 

F T V P f = V (det(F)p) - VF T : F" T det(F)p, (55) 

demonstrating the appearance of a hydrostatic stress-driven contribution to 3*. This is Darcy's Law for 
transport of a fluid down a pressure gradient. 

For the special case of a compressible, ideal fluid we have e f = e f (7y f ,p f ); i.e., the fluid stores strain 
energy as a function of its current, intrinsic density. Fluid saturation conditions hold in biological tissue, 
for which case the fluid volume fraction, / f , is simply the pore volume fraction. Recall from Section 
that the individual species deform with the common deformation gradient F. Therefore the pores deform 
homogeneously with the surrounding solid phase. Physically this corresponds to the pore size being smaller 
than the scale at which the homogenization assumption of a continuum theory holds. Momentarily ignoring 
changes in reference concentration of the fluid, we have F c = F. Then, since p^ = Po/ f , we can write 
e f (F, if , Pq) = e f (F, rf , p f f f ) = e f (rf , /5g/detF) = e f {rf , p f ). In this case a simple calculation shows that the 
hydrostatic pressure is 

_ p f de { 
P ~ ~dct(F) dp 1 ' 

and the stress divergence term is 

^' = -^'f) +VFT ^ v f' 

6.3 The Eshelby stress as a thermodynamic driving force 

Combining the stress divergence and chemical potential gradient contributions to the driving force for any 
species, and using the mass-specific Helmholtz free energy, ip L , we write, 

F T V -P'-Po (Ve' - 6Vrf) = V ■ (f t P 1 ) - VF T : P L - (56) 
Regrouping terms this expression is 

- V • (pW\el - F T P') + (Vp&) %\e - VF T : P\ (57) 
Eshelby stress, H 1 

Thus, the divergence of the well-known Eshelby stress tensor is also among the driving forces for mass 
transport. Also observe the presence of a strain gradient-dependent driving force, — VF T :P t in the devel- 
opments of Sections 16. 21 and 16.31 independent of the pressure gradient term for the fluid species. 
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Remark 1: The final version of the dissipation inequality (|50|l . and the mass balance equation can be 
manipulated to restrict the mathematic al form of the mass source. It is common to make the mass source 
depend upon the strain energy density I Harrigan and Hamilton! . 199.*^ while respecting the restriction 



im- 
posed by the dissipation inequality. This form is often used while modelling hard tissue. Such an approach 
leads to strain-mediated mass transport. However, with a strain-independent source, strain-mediated (or 
stress-mediated) mass transport would not be obtained with such a formulation. 

Remark 2: We expect that evaluation of the dissipation, D, using l|5U|l from field quantities in a bound- 
ary value problem will provide a test of soundness, and if necessary indications for improvement, of our 
constitutive models. 

Remark 3: Since soft biological tissues usually demonstrate rate-dependent response, it has been common 
to employ a solid viscoelastic constitutive model for them. This approach fits within our framework, with 
a modification of the internal energy to include its dependence upon internal variables that represent the 
viscoelastic stress-like parameters. However, a more physiologically-valid model may be one with a purely 
hyperelastic solid phase, and a viscous fluid. In such a composite model the rate-dependent behavior would 
arise from the fluid. 

Remark 4: The constitutive relations (|46|l and (|52() respectively specify the partial stress, P\ and flux, M L , 
of a species. The flux also implies the relative velocity, V L . The velocity of the solid phase, V is obtained 
from the local form of the balance of linear momentum for the system H21I) . With all these quantities known, 
the individual interaction forces between species, p^q 1 , can be obtained from (|26H . They are, however, not 
needed while solving for the balance of linear momentum of the system. 



7 A numerical example 

The theory developed in Sections EHEl has been implemented in a computational formulation, retaining 
much of the complexity of the coupled balance laws and constitutive relations. For realistic soft tissue 
material parameters, the contribution of the fluxes and interaction forces between species to the balance of 
linear momentum of the composite tissue is negligible. This simplification has been used. As a preliminary 
demonstration of the theory 7 , we present a computation of the coupled physics in the early stages of uniaxial 
extension of a cylindrical soft tissue specimen. The motivation for this model problem comes from our 
experimental model of engineered, functional tendon constructs grown in vitro, having the same cylindrical 
geometry. The exp erimental aspects of our broad-based project on soft tissue growth are described elsewhere 
1 Calve et all 12003^ . n addition to engineering scaffold-less tendon constructs from neonatal rat fibroblast 



cells, we have the ability to impose a range of mechanical, chemical, nutritional and electrical stimuli on 
them and study the tissue's response. Besides modelling these experiments, the mathematical formulation 
described here presents researchers with a vehicle for testing scenarios and framing hypotheses that can be 
experimentally-validated in our laboratory, thereby driving the experimental studies. 

7 This numerical section has been included mainly for completeness of this theoretical paper. A separate paper, currently in 
preparation, will present the computational formulation and contain a detailed examination of a number of initial and boundary 
value problems for growth. 
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7.1 Material models and parameters 



The engineered tendon construct is 12 mm in length and 1 mm 2 in area. In this paper an inter nal energy 
densit y for the solid phase based upo n the worm-like chain model is used. The reader is directed to 



(1997) and 



Bustamante et al 



Ricf ct al 



where the one-dimensional version of this model has been applied to long 
ch ain molecules. I t has been described and implemented into an anisotropic representative volume element 
by Bischoff et al. j2002^) . and is summarized here. The internal energy density of a single constituent chain 



of an eight-chain model (Figure is, 



Nk9 



r 

2L 



L 



4(1 -r/L) 



NkO ( 



4(1 - y/2A/L) 



Iog(A? A* X c 3 ) 



-2/3 



- l) + 2 7 l:£; c 



(58) 



Here, N is the density of chains, k is the Boltzmann constant, r is the end-to-end length of a chain, L is the 
fully-extended length, and A is the persistence length that measures the degree to which the chain departs 
from a straight line. The preferred orientation of tendon collagen is described by an anisotropic unit cell 
with sides a, b and c — see FigureOl All lengths in this model have been rendered non-dimensional (Table^) 
by dividing by the link length in a chain. 




Figure 3: Worm- like chains grouped into an initially anisotropic eight-chain model. 



The elastic stretches along the unit cell axes are respectively denoted Af , A| and A|, and E° = \{C C — 1) 
is the elastic Lagrange strain. The factors 7 and f3 control bulk compressibility. The end-to-end length is 
given by 



b 2 \f + c 2 Xf , X c j = y/Ni ■ C e Ni 



(59) 



Preliminary mechanical tests of the engineered tendon have been carried out in our laboratory but, at 
this stage, the worm-like chain model has not been calibrated to t hese tests. Instead, p ublished data for the 
worm-like chain, obtained by calibrating against rat cardiac tissue I Bischoff et al. . 2002]), has been employed. 
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The fluid phase was modelled as an ideal, nearly-incompressible fluid: 

p l e l {F c \pW) = i K (det(F° f ) - l) 2 , (60) 

where k is the fluid bulk modulus. 

Only a solid and a f luid p hase were included for the tissue. Low values were chosen for the mobilities of 



the fluid (Swa rtz et allll999ll with respect to the solid phase (see Tabled . In order to demonstrate growth, 
the solid phase must have a source term, II s (Section |2J, and the only other phase, the fluid, must have 
II f = —II s . Therefore, contrary to the case made in Section [3 a non-zero value of the fluid source, Il f , was 
assumed. A form motivated by first-order reactions was used: 

n f = -k l {pl -pIJ, n s = -n f , (6i) 

where k { is the reaction rate, and p^. . is the initial fluid concentration. This term acts as a source for the 
solid when p { > pr Q _ . , and a sink when p Q < p . 

In a very simple approximation, the fluid's mixing entropy was written as 



Recall that in the notation of Section [21 M f is the fluid's molecular weight. 



(62) 



Table 1: Material parameters used in the analysis 



Parameter 


Symbol 


Value 


Units 


Chain density 


N 


7 x 10 21 


_3 

m 


Temperature 





310.0 


K 


Persistence length 


A 


1.3775 




Fully-stretched length 


L 


25.277 




Unit cell axes 


a, 6, , c 


9.2981, 12.398, 6.1968 




Bulk compressibility factors 


7, 


1000, 4.5 




Fluid bulk modulus 


K 


1 


GPa 


Fluid mobility tensor components 


Dn, D22, -D33 


1 x 10~ 8 , 1 x 10~ 8 , 1 x 10~ 8 


m~ 2 sec 


Fluid conversion reaction rate 


fc f 


-1. x 10~ 7 


sec -1 


Gravitational acceleration 


9 


9.81 


m.sec -2 


Molecular weight of fluid 


M f 


2.9885 x 10" 23 


kg 



7.2 Boundary and initial conditions; coupled solution method 

Boundary conditions for mass transport consisted of the specified fluid concentration at all external surfaces 
of the cylinder. This value was fixed at 500 kg. m -3 . With these boundary conditions the fluid flux normal 
to surfaces of the specimen is determined by solving the initial and boundary value problem. The bottom 
planar surface was fixed in the e% direction and a displacement was applied at the top surface, also in the e% 
direction, to give a nominal strain rate of 0.05 sec -1 in the e 3 direction. This is the only mechanical load on 
the problem. Initial conditions were p { (X,0) = 500 kg. m -3 , Pq(X,0) — 500kg.nU 3 , and for the mechanical 
problem, u(X,0) = 0, V(X,0) = 0. 
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The coupled probl em was solved by a staggered scheme based upon operator splits I Armerol 



Garikipati et all 12001). The details will be presented in a future communication that will focus upon com- 
putational aspects and numerical examples. Here we only mention that the staggered scheme consists of 
identifying the displacement and species concentrations as primitive variables associated with the mechani- 
cal and mass transport problems. The mechanical problem is solved holding the concentrations fixed. The 
resulting displacement field is then held constant to solve the mass transport problem. The transient so- 
lution is obtain ed for mechanics using energy-momentum conserving schemes I Simo and Tarnow , 1992albl : 
Gonzalezlll99^ . and for mass transport using the B ackward Euler Me thod. Hexahedral elements are em- 
ployed, combined with nonlinear projection methods Ipimo et all Il985|) to treat the near-incompressibility 
imposed by the fl uid. The num erical formulation has been implemented within the nonlinear finite element 
program, FEAP |Tavlorl fl999^ . 



7.3 Evolution of stress and flux 

The following contour plots represent the stress, and various contributions to the total flux in the early stages 
of loading of the model problem. Symmetry has been employed to model a single quadrant of the cylinder. 

The longitudinal stress, 033 in Figure 0] arises from the stretch and the evolution in concentration. 




I 



7.61 E-01 
7.50E+01 
6.00E+01 
4.50E+01 
3.00E+01 
1 .50E+01 
0.00E+OO 
-1.50E+01 
-3.00E+01 
-4.50E+01 
-6.00E+01 
-7.50E+01 
-7.46E-02 



Time= 1.00E-09 




1 



7.63E+01 
7.50E+01 
6.00E+01 
4.50E+01 
3.00E+01 
1.50E+01 
0.00E+00 
-1 .50E+01 
-3.00E+01 
-4.50E+01 
-6.00E+01 
-7.50E+01 
-9.63E+00 



Time = 1 .OOE-07 



(a) (b) 

Figure 4: Longitudinal Cauchy stress, 033 (Pa) at 1 nanosec. and 100 nanosec. after the beginning of loading. 



The flux contributions in Figures 151- 1101 can be summarized as follows: The fluid flux is dominated by 
the contribution from the stress gradient in the direction. The latter arises as the stress (033) wave of 
tension travels down the cylinder in the first few microseconds after application of the load (the time taken 
to travel the length of the cylinder is 12/xsec). Additionally, as the fluid concentration changes due to the 
flux, it causes a further change in the stress (Section^. Other flux terms are qualitatively sensible; i.e., 
their directions are consistent with the physics of the problem, as argued in each of the figure captions 8 . 
There is some loss of axial symmetry in a few of the plots due to the coarseness of the finite element mesh 
for this example. It appears that spatial oscillations in the solution lead to a further loss of symmetry in 

8 In order to compare the flux contributions, they have all been plotted on the same scale: —1 X 10~ 4 1 X kg.m - 2 sec. 

However the plots also show the maximum and minimum field values at the top and bottom of the legend bars. These values 
represent a better comparison of the relative flux magnitudes. 
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Figure 5: Stress gradient-driven flux (kg.m _2 sec) in the e% direction at 1 nanosec. and lOOnanosec. after 
the beginning of loading. The positive values indicate an upward flux corresponding to a tensile 033 wave 
travelling downwards. 



Figures [fj>-[7|3 and If Ufa,. These oscillations arise due to large and dominant advective terms, and need to be 
remedied by stabilized finite element methods. Here, we only aim to demonstrate that various driving forces 
for mass transport are in agreement with their theoretical underpinnings in the paper. The resorption of the 
solid phase is shown indirectly in Figure ITT1 A positive fluid source, II , means that II s < 0. Since II s is the 
only term balancing dp^/dt [see J2J], it follows that dp^/dt < 0. 



8 Discussion and conclusions 



A general framework for growth of biological tissue has been presented in this paper. While simplified 
models were used for source terms, (II 1 ), they can be formulated on the basis of the kinetics of chemical 
reactions in order to develop more realistic growth laws. This approach, we believe, is fundamental to 
a proper treatment of mass transport in tissue. Results are obtained that differ fundamentally from the 
classical setting of continuum mechanics. Most notable among these differences are the mass fluxes driven 
by gradients in stress, strain, energy and entropy, in addition to body force and inertia. Importantly, though 
our treatment differs from classical mixture theory, the two are fully consistent as established at several points 
in this paper. The balance laws in Section [2] and [3] introduce a degree of coupling between the phenomena 
of mass transport and mechanics. This is visible most transparently in the balance of linear momentum (|26|) 
that includes mass fluxes, M L . The balance of mass, described by Equations J2J) and l@, is also dependent 
upon the mechanics as the discussion in Section l6~Tl makes clear. Notably, this ensures mechanics- mediated 
mass transport even with a mass source that is independent of strain/stress, as the discussion at the end 
of Section 16.31 establishes. The discussion in Sections I6.IH6.3I provides many insights into the nature of 
this coupling. The mechanics problem also has a constitutive dependence upon mass concentration, via 
(14611 and the fact that the growth deformation gradient tensor, F s is determined by the concentration. The 
viscoelastic nature of the composite tissue would emerge naturally from a model incorporating a hyperelastic 
solid and viscous fluid. 

We have formally allowed all species to be load bearing and develop a stress. At the scales that are 
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Figure 6: Internal energy gradient-driven flux, (kg.m _2 sec) in the e% direction at 1 nanosec. and 100 nanosec. 
after the beginning of loading. The positive values indicate an upward flux. This corresponds to a lower 
energy near the top of the cylinder as the tensile stress (0-33) wave travels downward and relaxes some of the 
strain energy of contraction. 



of interest in a tissue, the only relevant load bearing species are the solid and fluid phases. Nevertheless, 
inasmuch as a transported species such as a nutrient has a molecular structure that can be subject to loads 
at the scale of pico-newtons, it is not inconsistent to speak of the partial stress of this species. Since the 
constitutive relation H46|) indicates that the partial stresses are scaled by concentrations, the contribution to 
total stress from any species besides the solid and fluid phases will be negligible. 

We have chosen to leave remodelling out of our formulation in this communication, to focus upon the 
above issues. Remodelling includes any evolution in properties, state of stress, material symmetry, volume 
or shape brought about by microstructural changes. In the development of biological tissue, growth and 
remodelling occur simultaneously. As density changes due to growth, the material also remodels by mi- 
crostructural evolution within the neighborhood of each point. A rigorous tre atment of this phenom enon 



has been presented in the continuum mechanical setting in a companion paper ijGarikipati et al 
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Figure 9: Concentration gradient-driven flux (kg.m _2 sec) in the e% direction at 1 nanosec. and 100 nanosec. 
after the beginning of loading. Note that the maximum and minimum values are many orders of magnitude 
lower than for the other flux contributions reported above. This is a demonstration of mechanics influences 
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Figure 10: Total flux (kg.m 2 sec) in the direction at 1 nanosec. and lOOnanosec. after the beginning of 
loading. The positive values indicate an upward flux, dominated by the stress gradient driven contribution. 
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Figure 11: Rate of fluid production, II f (kg.m 3 .sec 1 ), at 1 nanosec. and lOOnanosec. after the beginning of 
loading. The positive values indicate that the local fluid concentrations have fallen below their initial values. 
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